! idx class for QuickPIC Open Source 1.0
! update: 04/18/2016

module idx_class
use parallel_pipe_class
use input_class


implicit none

private

public :: idx

type idx


class(parallel_pipe), pointer, public :: p => null()

!glb_vtx global meshgrid range, coarse layer reference frame
!glb_id global meshgrid real id range, glb_id * dx
!loc_id local meshgrid real id range
!fld_id field id that combines the true field
! rdx rdy rdz are relative cell size relative to the elementary coarse level(not the coarsest in solver)
! dx dy dz, are real cell size .

! parallel parameter
! noff = smallest global gridpoint in y, noff + loc_id_y = glb_id_y
! nx, ny size of global array data in x and y direction
! nxp, nyp size of local array data in x and y direction

integer, dimension(2) :: glb_vtx_x , glb_vtx_y
integer, dimension(2) :: glb_id_x , glb_id_y
integer, dimension(2) :: loc_id_x, loc_id_y
integer, dimension(2) ::  fld_id_x, fld_id_y
integer, dimension(2) :: push_id_x, push_id_y
integer, dimension(2) :: glb_id_z, loc_id_z
!this is added for we want to push particle at edge
integer, dimension(2) :: dep_loc_x, dep_loc_y, dep_glb_x, dep_glb_y, dep_glb_z

integer, dimension(2) :: stage_id_z, fld_id_z, dep_loc_z
real, dimension(2) :: save_vtx_x, save_vtx_y,save_vtx_z
real, dimension(2):: box_x, box_y, box_z
real, dimension(2) :: dep_box_x, dep_box_y, dep_box_z
real, dimension(2) :: push_box_x, push_box_y
logical, dimension(4) :: is_phys_bnd
logical, dimension(2) :: is_phys_bnd_z = (/ .false., .false. /)
integer :: grid_level, sublevel
!level index and subcycle
integer, dimension(2):: grid_multi

real :: rdx, rdy, rdz
real :: dx, dy, dz
real :: posoffx, posoffy, posoffz
!parallel group number
! integer :: pg


!!!!!!!!!!!!!!!!!!!!Need to change level into logical
!level id
integer :: lid,level0
integer :: level

integer :: noff, noffz
integer :: glb_nx, glb_ny, glb_nz, loc_nx, loc_ny, fld_nx, fld_ny,loc_nz,fld_nz
integer :: nproc_y

character(len=12) :: fld_name
character(len=15) :: fpart_name
!p1_l1_s1 means processor 1 level 1 subcycle 1


!compatible variable
  integer :: state, gcells

contains
    generic :: id2pos_x => id2pos_x1, id2pos_x2
    generic :: id2pos_y => id2pos_y1, id2pos_y2
    generic :: id2pos_z => id2pos_z1, id2pos_z2
    generic :: gen => gen_dependent_parameter
    generic :: cover =>cover2d, cover3d, cover2d_i
    generic :: dep_cover =>dep_cover2d, dep_cover3d
    generic :: push_cover => push_cover2d
    ! generic :: glb2loc => glb2loci, glb2locr


    procedure :: print_idx,print_idx_compact
    procedure :: id_loc2glb, id_glb2loc,gen_dependent_parameter,id_loc2glb_z
    procedure :: id2pos_x1, id2pos_y1,id2pos_x2, id2pos_y2
    procedure :: id2pos_z1, id2pos_z2
    procedure :: pos2id_x, pos2id_y, pos2id_z
    procedure :: inbox, cover3d,cover2d, cover2d_i, cover2d_in
    procedure :: dep_cover2d, dep_cover3d, push_cover2d
    procedure :: copy_idx
    procedure :: coverz

end type


contains
!

subroutine print_idx(this)

  implicit none

  class(idx), intent(inout) :: this

  print *,this%fpart_name,"noff", this%noff
  print *,this%fpart_name,"noffz", this%noffz
  print *,this%fpart_name,"grid_level", this%grid_level
  print *,this%fpart_name,"glb_id_x", this%glb_id_x
  print *,this%fpart_name,"glb_id_y", this%glb_id_y
  print *,this%fpart_name,"glb_vtx_x", this%glb_vtx_x
  print *,this%fpart_name,"glb_vtx_y", this%glb_vtx_y
  print *,this%fpart_name,"save_vtx_x", this%save_vtx_x
  print *,this%fpart_name,"save_vtx_z", this%save_vtx_z

  print *,this%fpart_name,"loc_id_x", this%loc_id_x
  print *,this%fpart_name,"loc_id_y", this%loc_id_y
  print *,this%fpart_name,"dep_loc_x", this%dep_loc_x
  print *,this%fpart_name,"dep_loc_y", this%dep_loc_y
  print *,this%fpart_name,"dep_glb_x", this%dep_glb_x
  print *,this%fpart_name,"dep_glb_y", this%dep_glb_y
  print *,this%fpart_name,"dep_glb_z", this%dep_glb_z
  print *,this%fpart_name,"box_x", this%box_x
  print *,this%fpart_name,"box_y", this%box_y
  print *,this%fpart_name,"box_z", this%box_z
  print *,this%fpart_name,"dep_box_x", this%dep_box_x
  print *,this%fpart_name,"dep_box_y", this%dep_box_y
  print *,this%fpart_name,"dep_box_z", this%dep_box_z
  print *,this%fpart_name,"fld_id_y", this%fld_id_y
  print *,this%fpart_name,"fld_id_x", this%fld_id_x
  print *,this%fpart_name,"glb_id_z", this%glb_id_z
  print *,this%fpart_name,"grid_multi", this%grid_multi
  print *,this%fpart_name,"level", this%level
  print *,this%fpart_name,"rdx", this%rdx
  print *,this%fpart_name,"rdy", this%rdy
  print *,this%fpart_name,"rdz", this%rdz

  print *,this%fpart_name,"glb_nx", this%glb_nx
  print *,this%fpart_name,"glb_ny", this%glb_ny
  print *,this%fpart_name,"loc_nx", this%loc_nx
  print *,this%fpart_name,"loc_ny", this%loc_ny
  print *,this%fpart_name,"fld_nx", this%fld_nx
  print *,this%fpart_name,"fld_ny", this%fld_ny

  print *,this%fpart_name,"is_phys_bnd", this%is_phys_bnd

end subroutine print_idx

subroutine print_idx_compact(this)

  implicit none

  class(idx), intent(inout) :: this

  print *,this%fpart_name,"noff", this%noff
  print *,this%fpart_name,"grid_level", this%grid_level
  print *,this%fpart_name,"glb_id_x", this%glb_id_x
  print *,this%fpart_name,"glb_id_y", this%glb_id_y
  print *,this%fpart_name,"loc_id_x", this%loc_id_x
  print *,this%fpart_name,"loc_id_y", this%loc_id_y
  print *,this%fpart_name,"dep_loc_x", this%dep_loc_x
  print *,this%fpart_name,"dep_loc_y", this%dep_loc_y
  print *,this%fpart_name,"dep_loc_z", this%dep_loc_z
  print *,this%fpart_name,"dep_glb_z", this%dep_loc_z
  print *,this%fpart_name,"stage_id_z", this%stage_id_z
  print *,this%fpart_name,"fld_id_z", this%fld_id_z
  print *,this%fpart_name,"is_phys_bnd_z", this%is_phys_bnd_z
  print *,this%fpart_name,"noffz", this%noffz

  print *,this%fpart_name,"fld_id_y", this%fld_id_y
  print *,this%fpart_name,"fld_id_x", this%fld_id_x
  print *,this%fpart_name,"glb_id_z", this%glb_id_z
  print *,this%fpart_name,"grid_multi", this%grid_multi
  print *,this%fpart_name,"level", this%level
  print *,this%fpart_name,"rdx", this%rdx
  print *,this%fpart_name,"rdy", this%rdy
  print *,this%fpart_name,"rdz", this%rdz

end subroutine print_idx_compact

subroutine copy_idx(this,id)

  implicit none

  class(idx), intent(inout) :: this
  class(idx), intent(in) :: id

  this%grid_level = id%grid_level
  ! this%pg = id%pg
  this%sublevel = id%sublevel
  this%p => id%p
  this%level = id%level
  this%glb_vtx_x = id%glb_vtx_x
  this%glb_vtx_y = id%glb_vtx_y
  this%grid_multi = id%grid_multi
  this%glb_id_x = id%glb_id_x
  this%glb_id_y = id%glb_id_y
  this%glb_id_z = id%glb_id_z

  this%save_vtx_x = id%save_vtx_x
  this%save_vtx_y = id%save_vtx_y
  this%save_vtx_z = id%save_vtx_z

  this%lid = id%lid
  this%level0 = id%level0
  this%rdx = id%rdx
  this%rdy = id%rdy
  this%rdz = id%rdz

  this%dx = id%dx
  this%dy = id%dy
  this%dz = id%dz

  this%glb_nx = id%glb_nx
  this%glb_ny = id%glb_ny
  this%glb_nz = id%glb_nz


  this%nproc_y = id%nproc_y
  this%noff = id%noff
  this%noffz = id%noffz
  this%loc_id_x = id%loc_id_x
  this%loc_id_y = id%loc_id_y
  this%fld_id_y = id%fld_id_y
  this%fld_id_x = id%fld_id_x
  this%fld_id_z = id%fld_id_z
  this%loc_nx = id%loc_nx
  this%loc_ny = id%loc_ny
  this%fld_ny = id%fld_ny
  this%fld_nx = id%fld_nx
  this%fld_nz = id%fld_nz

  this%dep_glb_x = id%dep_glb_x
  this%dep_glb_y = id%dep_glb_y
  this%dep_glb_z = id%dep_glb_z
  this%dep_loc_x = id%dep_loc_x
  this%dep_loc_y = id%dep_loc_y
  this%dep_loc_z = id%dep_loc_z

  this%posoffx = id%posoffx
  this%posoffy = id%posoffy
  this%posoffz = id%posoffz

  this%box_x = id%box_x
  this%box_y = id%box_y
  this%box_z = id%box_z
  this%dep_box_x = id%dep_box_x
  this%dep_box_y = id%dep_box_y
  this%dep_box_z = id%dep_box_z

  this%push_id_x = id%push_id_x
  this%push_id_y = id%push_id_y

  this%push_box_x = id%push_box_x
  this%push_box_y = id%push_box_y

  this%is_phys_bnd = id%is_phys_bnd
  this%is_phys_bnd_z = id%is_phys_bnd_z

  this%stage_id_z = id%stage_id_z

end subroutine copy_idx

subroutine gen_dependent_parameter(this, bd_width)

  implicit none

  class(idx), intent(inout) :: this
  integer, intent(in) :: bd_width
  integer :: i

    this%glb_nx = this%glb_id_x(2) - this%glb_id_x(1) + 1
    this%glb_ny = this%glb_id_y(2) - this%glb_id_y(1) + 1
    this%glb_nz = this%glb_id_z(2) - this%glb_id_z(1) + 1
    this%loc_nx = this%loc_id_x(2) - this%loc_id_x(1) + 1
    this%loc_ny = this%loc_id_y(2) - this%loc_id_y(1) + 1
    this%fld_nx = this%fld_id_x(2) - this%fld_id_x(1) + 1
    this%fld_ny = this%fld_id_y(2) - this%fld_id_y(1) + 1
    this%fld_nz = this%fld_id_z(2) - this%fld_id_z(1) + 1




    if (this%grid_level == 1) then

      this%push_id_x(1) = this%glb_id_x(1)
      this%push_id_x(2) = this%glb_id_x(2) + 1


      ! this%push_id_y(1) = this%glb_id_y(1)
      ! this%push_id_y(2) = this%glb_id_y(2)

      this%push_id_y(1) = this%glb_id_y(1)
      this%push_id_y(2) = this%glb_id_y(2)

    else
      this%push_id_x(1) = this%glb_id_x(1) + bd_width
      this%push_id_x(2) = this%glb_id_x(2) - bd_width
      this%push_id_y(1) = this%glb_id_y(1) + bd_width
      this%push_id_y(2) = this%glb_id_y(2) - bd_width


      ! this%push_id_x(1) = this%glb_id_x(1)
      ! this%push_id_x(2) = this%glb_id_x(2)
      ! ! this%push_id_y(1) = this%glb_id_y(1)
      ! ! this%push_id_y(2) = this%glb_id_y(2)
      ! this%push_id_y(1) = this%glb_id_y(1)
      ! this%push_id_y(2) = this%glb_id_y(2)
    end if

    ! if (this%grid_level == 1) then
    !   this%push_id_x(2) = this%push_id_x(2) + 1
    !   this%push_id_y(2) = this%push_id_y(2) + 1
    ! end if


    do i = 1, 2
      !!!!!!!!!!!!!!!!!!!

      this%box_x(i) = this%id2pos_x(this%glb_id_x(i))
      this%box_y(i) = this%id2pos_y(this%glb_id_y(i) - this%noff)
      this%box_z(i) = this%id2pos_z(this%glb_id_z(i))
      this%dep_box_x(i) = this%id2pos_x(this%dep_glb_x(i))
      this%dep_box_y(i) = this%id2pos_y(this%dep_glb_y(i) - this%noff)
      this%dep_box_z(i) = this%id2pos_z(this%dep_glb_z(i))
      this%push_box_x(i) = this%id2pos_x(this%push_id_x(i))
      this%push_box_y(i) = this%id2pos_y(this%push_id_y(i) - this%noff)

    end do
    if (this%grid_level == 1) then


      this%box_x(2) = this%box_x(2) + 1.0
      this%box_y(2) = this%box_y(2) + 1.0
      this%box_z(2) = this%box_z(2) + 1.0

      this%dep_box_x(1) = this%dep_box_x(1) - 1.0
      this%dep_box_y(1) = this%dep_box_y(1) - 1.0
      this%dep_box_z(1) = this%dep_box_z(1) - 1.0
      ! this%push_id_x(1) = this%push_id_x(1) - 1
      ! this%push_id_y(1) = this%push_id_y(1) - 1

    end if

end subroutine gen_dependent_parameter

function id_loc2glb(this, loc_id)
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: loc_id
  integer :: id_loc2glb

  id_loc2glb = loc_id + this%noff
end function id_loc2glb

function id_loc2glb_z(this, loc_id_z)
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: loc_id_z
  integer :: id_loc2glb_z

  id_loc2glb_z = loc_id_z + this%noffz
end function id_loc2glb_z

function id_glb2loc(this, glb_id)
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: glb_id
  integer :: id_glb2loc

  id_glb2loc = glb_id - this%noff
end function id_glb2loc


function pos2id_x(this, pos_x)
  implicit none
  class(idx), intent(in) :: this
  real, intent(in):: pos_x
  real :: pos2id_x

  pos2id_x = (pos_x - this%posoffx)/this%rdx

end function pos2id_x

function pos2id_y(this, pos_y)
  implicit none
  class(idx), intent(in) :: this
  real, intent(in):: pos_y
  real :: pos2id_y

  pos2id_y = (pos_y - this%posoffy)/this%rdy - this%noff

end function pos2id_y

function pos2id_z(this, pos_z)
  implicit none
  class(idx), intent(in) :: this
  real, intent(in):: pos_z
  real :: pos2id_z

  pos2id_z = (pos_z - this%posoffz)/this%rdz - this%noffz

 ! print *, pos2id_z, pos_z, this%posoffz, this%rdz

end function pos2id_z

function id2pos_x1(this, id_x)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: id_x
  real :: id2pos_x1

  id2pos_x1 = id_x * this%rdx + this%posoffx

end function id2pos_x1

function id2pos_y1(this, id_y)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: id_y
  real :: id2pos_y1

  id2pos_y1 = (id_y + this%noff) * this%rdy + this%posoffy

end function id2pos_y1

function id2pos_x2(this, id_x)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: id_x
  real :: id2pos_x2

  id2pos_x2 = id_x * this%rdx + this%posoffz

end function id2pos_x2

function id2pos_y2(this, id_y)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: id_y
  real :: id2pos_y2

  id2pos_y2 = (id_y + this%noff) * this%rdy + this%posoffy

end function id2pos_y2

function id2pos_z1(this, id_z)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: id_z
  real :: id2pos_z1

  id2pos_z1 = id_z * this%rdz + this%posoffz

end function id2pos_z1

function id2pos_z2(this, id_z)
  ! transfer loc id(can be real number) to global position assuming elementary level d = 1
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: id_z
  real :: id2pos_z2

  id2pos_z2 = id_z * this%rdz + this%posoffz

end function id2pos_z2

function inbox(this,lpos)
  !return if the lpos(the \xi index) is in this level.
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: lpos
  logical :: inbox
  if (this%level == 0) then
    inbox = .false.
  else if ((lpos + this%noffz) .ge. this%glb_id_z(1) .and. (lpos+ this%noffz) .le. this%glb_id_z(2)) then
    inbox = .true.
  else
    inbox = .false.
  end if

end function inbox

function cover3d(this,x,y,z)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: x,y,z
  logical :: cover3d

  if (this%level == 0 ) then
    cover3d = .false.
  else if (x .ge. this%box_x(1) .and. x .le. this%box_x(2) .and.&
  &  y .ge. this%box_y(1) .and. y .le. this%box_y(2) .and. &
  &  z .ge. this%box_z(1) .and. z .le. this%box_z(2) ) then
    cover3d= .true.
  else
    cover3d = .false.
  end if

end function cover3d


function cover2d(this,x,y)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: x,y
  logical :: cover2d


  if (this%level == 0) then
    cover2d = .false.
  else if (x .ge. this%box_x(1) .and. x .le. this%box_x(2) .and.&
  &  y .ge. this%box_y(1) .and. y .le. this%box_y(2)) then
    cover2d= .true.
  else
    cover2d = .false.
  end if

end function cover2d

function cover2d_i(this,x,y)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: x,y
  logical :: cover2d_i


  if (this%level == 0) then
    cover2d_i = .false.
 else if (x .ge. this%glb_id_x(1) .and. x .le. this%glb_id_x(2) .and.&
  &  this%id_loc2glb(y) .ge. this%glb_id_y(1) &
  & .and. this%id_loc2glb(y) .le. this%glb_id_y(2)) then
    cover2d_i = .true.
  else
    cover2d_i = .false.
  end if

end function cover2d_i

function cover2d_in(this,x,y)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: x,y
  logical :: cover2d_in


  if (this%level == 0) then
    cover2d_in = .false.
  else if (x > this%glb_id_x(1) .and. x < this%glb_id_x(2) .and.&
  &  this%id_loc2glb(y) > this%glb_id_y(1) &
  & .and. this%id_loc2glb(y) < this%glb_id_y(2)) then

    cover2d_in = .true.
  else
    cover2d_in = .false.
  end if

end function cover2d_in

function dep_cover3d(this,x,y,z)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: x,y,z
  logical :: dep_cover3d

  if (this%level == 0 ) then
    dep_cover3d = .false.
  else if (x .ge. this%dep_box_x(1) .and. x .le. this%dep_box_x(2) .and.&
  &  y .ge. this%dep_box_y(1) .and. y .le. this%dep_box_y(2) .and. &
  &  z .ge. this%dep_box_z(1) .and. z .le. this%dep_box_z(2) ) then
    dep_cover3d= .true.
  else
    dep_cover3d = .false.
  end if

end function dep_cover3d


function dep_cover2d(this,x,y)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: x,y
  logical :: dep_cover2d


  if (this%level == 0) then
    dep_cover2d = .false.
  else if (x .ge. this%dep_box_x(1) .and. x .le. this%dep_box_x(2) .and.&
  &  y .ge. this%dep_box_y(1) .and. y .le. this%dep_box_y(2)) then
    dep_cover2d= .true.
  else
    dep_cover2d = .false.
  end if

end function dep_cover2d

function push_cover2d(this,x,y)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  real, intent(in) :: x,y
  logical :: push_cover2d


  if (this%level == 0) then
    push_cover2d = .false.
  else if (x .ge. this%push_box_x(1) .and. x .le. this%push_box_x(2) .and.&
  &  y .ge. this%push_box_y(1) .and. y .le. this%push_box_y(2)) then
    push_cover2d= .true.
  else
    push_cover2d = .false.
  end if

end function push_cover2d

function coverz(this,pos_z)
  !return if particle at position x,y,z is covered in this level of mesh, the unit is normalized by level_0 dx,dy,dz
  implicit none
  class(idx), intent(in) :: this
  integer, intent(in) :: pos_z
  logical :: coverz


  if (((pos_z  + this%noffz) > this%glb_id_z(2) &
  & .or. (pos_z  + this%noffz) < this%glb_id_z(1)) .and.  pos_z > 0 ) then
    coverz = .false.
  else
    coverz = .true.
  end if

end function coverz

end module idx_class
